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Abstract. In an effort to study the applicability of adaptive mesh refinement (AMR) techniques 
to atmospheric models an interpolation-based spectral element shallow water model on a cubed- 
sphere grid is compared to a block-structured finite volume method in latitude-longitude geometry. 
Both models utilize a non-conforming adaptation approach which doubles the resolution at fine-coarse 
mesh interfaces. The underlying AMR libraries are quad-tree based and ensure that neighboring 
regions can only differ by one refinement level. The models are compared via selected test cases from 
a standard test suite for the shallow water equations. They include the advection of a cosine bell, a 
steady-state geostrophic flow, a flow over an idealized mountain and a Rossby-Haurwitz wave. Both 
static and dynamics adaptations are evaluated which reveal the strengths and weaknesses of the 
AMR techniques. Overall, the AMR simulations show that both models successfully place static and 
dynamic adaptations in local regions without requiring a flne grid in the global domain. The adaptive 
grids reliably track features of interests without visible distortions or noise at mesh interfaces. Simple 
threshold adaptation criteria for the geopotential height and the relative vorticity are assessed. 

1. Introduction. Simulating the climate is a grand challenge problem requiring 
multiple, century-long integrations of the equations modeling the Earth's atmosphere. 
Consequently, grid resolutions in atmospheric climate models are much coarser than 
in numerical weather models, where accurate predictions are limited to about ten 
days. However, it has been recognized that localized flow structures, like hurricanes, 
may play an important role in obtaining the correct climate signal. In particular, 
there exists a requirement for localized mesh refinement in regional climate modeling 
studies. It is also conjectured, at this time, that the lack of resolution in the tropics is 
the cause of the inability of most climate models to capture the statistical contribution 
of extreme weather events. 

While mesh adaptation is a mature field in computational fluid dynamics, cur- 
rently, the only fully operational adaptive weather and dispersion model is OMEGA 
[U [21] . The latter is based on the finite- volume approach and uses conforming Delau- 
nay meshes that are locally modified and smoothed. Researchers involved in hurricane 
or typhoon predictions where amongst the very first to experiment with variable reso- 
lutions in atmospheric models. The preferred method consisted of nesting finer meshes 
statically into coarser ones [331132] . Nesting is a wide spread technique in current mod- 
els, especially mesoscale models, to achieve high resolutions [14] in local regions. Truly 
adaptive models were first developed by [5T], [SD] and [TU]. [SD] based their adaptive 
mesh refinement (AMR) strategy on a truncation error estimate [49| . Continuing in 
the realm of second order conforming methods [24] developed a dynamically adaptive 
Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) [25] . 
MPDATA was thoroughly reviewed in [52] because of its various numerical qualities. 
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The discussion addressed generalized space-time coordinates that enable continuous 
mesh deformations [IS] as well as the generalization of MPDATA to unstructured 
grids [S3] . Recently, adaptive schemes on the sphere were also discussed by [IS] , [51] , 
[3] and |36| where only |26j utilizes a non-conforming approach. A review of adaptive 
methods in atmospheric modeling is given in [2] . 

It is recognized that high-order spatial resolution is necessary in climate modeling. 
To paraphrase [5], spectral methods like continuous Galerkin (CG) or discontinuous 
Galerkin (DG) are blessedly free of the spurious wave dispersion induced by low-order 
methodf0. The corresponding reduction in grid points directly lowers the amount of 
costly column physics evaluations. Multiple efforts to include conforming mesh re- 
finement techniques in models with high-order numerical methods were initiated |17| , 
[l3] . [20] . [19] and [47]. However, supporting conforming adaptive meshes has nega- 
tive effects on the Courant-Friedrich-Levy (CFL) stability condition. One solution is 
to consider non-conforming elements or blocks. 

In this paper two adaptive mesh techniques for 2D shallow water flows on the 
sphere are compared. In particular, the study focuses on the characteristics of the 
AMR approach in the cubed-sphere spectral element model (SEM) by [54] and com- 
pares it to the adaptive finite volume (FV) method by ^H]- Both SEM and FV are 
AMR models of the non-conforming type. A thorough study of the dynamically adap- 
tive FV scheme based on quadrilateral control volumes can also be found in [25] and 
[27j . In particular, 2D shallow water and 3D primitive-equation numerical experi- 
ments were conducted using a latitude-longitude grid on the sphere. This adaptive 
model is built upon the finite volume technique by [37] and [38]. The spectral ele- 
ment model is based on the ideas of [H]. The non-conforming treatment follows the 
interpolation procedure by [12j . The spectral element method was originally devel- 
oped for incompressible fluid flows. Meanwhile, it has been adapted by many authors 
[22l [57l [m [58] for global atmospheric and oceanic simulations. 

The paper is organized as follows. In Section [2] the shallow water equations 
are introduced which are the underlying equation for our model comparison. Both 
models SEM and FV are briefly reviewed in Section [2] This includes a discussion 
of the adaptive mesh approach for the spectral elements in SEM and the latitude- 
longitude blocks in FV. In Section U] the characteristics of the AMR techniques are 
tested using selected test cases from the [63] shallow water test suite. They include 
the advection of a cosine bell, a steady-state geostrophic flow, a flow over a mountain 
and a Rossby-Haurwitz wave. The findings are summarized in Section [5] 

2. Shallow water equations. The shallow water equations have been used as 
a test bed for promising numerical methods by the atmospheric modeling community 
for many years. They contain the essential wave propagation mechanisms found in 
atmospheric General Circulation Models (GCMs). The linearized primitive equations 
yield a series of layered shallow water problems where the mean depth of each layer 
is related to the maximum wave speed supported by the medium. These are the 
fast-moving gravity waves and nonlinear Rossby waves. The latter are important 
for correctly capturing nonlinear atmospheric dynamics. The governing equations of 
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motion for the inviscid flow of a free surface are given by 

1^ + (/ + k X V + iv (vv) + V($ + $.) = (2.1) 

— +V.($v) = (2.2) 

where t is the time, $ = gh symbohzes the geopotential, h is the height above sea level 
and g denotes the gravitational acceleration. Furthermore, <I>s = ghg is the surface 
geopotential, hs symbolizes the height of the orography, v stands for the horizontal 
velocity vector, k denotes the radial outward unit vector, and / and C are the Coriolis 
parameter and relative vorticity, respectively. V indicates the horizontal gradient 
operator, whereas V- stands for the horizontal divergence operator. The momentum 
equation (j2.ip is written in its vector-invariant form. 

3. Description of the adaptive shallov^f v^fater models. The spectral ele- 
ment method (SEM) is a combination of ideas coming from the finite-element method 
and of the spectral method. The order of accuracy is determined by the degree 
of the local basis functions within a finite element. The basis consists of Lagrange 
polynomials passing through Gauss-Legendre-Lobatto (GLL) or Gauss-Legendre (GL) 
quadrature points facilitating greatly the evaluation of the integrals appearing in a 
weak formulation. For the tests here, fifth degree basis functions for scalars are em- 
ployed which utilize 6x6 GL quadrature points per spectral element. The SEM grid is 
based on a projection of a cube inscribed into a sphere ([IH]), a so-called cubed-sphere 
grid which consists of 6 faces. These are further subdivided into spectral elements 
quadrangles which are evenly distributed over the surface of the sphere in the non- 
adapted case. Note that the numerical approach in SEM is non-monotonic and does 
not conserve mass. Nevertheless, the variation in the total mass is small over typical 
forecast periods of two weeks. 

This is in contrast to the monotonic and mass-conservative FV model which was 
originally developed by [38]. It utilizes the third order Piecewise Parabolic Method 
(PPM) [7] which was first designed for compressible fluids with strong shocks. The 
PPM algorithm applies monotonicity constraints that act as a nonlinear scale-selective 
dissipation mechanism. This dissipation primarily targets the flow features at the 
smallest scales. 

Both shallow water models arc characterized in more detail below. In particular, 
the individual adaptive mesh approaches are described which implement adaptations 
in the horizontal directions. The time step, on the other hand, is not adapted, except 
for selected advection tests. Therefore, the chosen time step must be numerically 
stable on the finest grid in an adapted model run. 

3.1. Spectral element (SEM) shallow water model. 

3.1.1. Curvilinear coordinates: cubed sphere. The flux form shallow- water 
equations in curvilinear coordinates arc described in [48j . Let ai and a.2 be the 
covariant base vectors of the transformation between inscribed cube and spherical 
surface. The metric tensor of the transformation is dcflned as Gij = SLi ■ a.j. Covariant 
and contravariant vectors are related through the metric tensor by Ui = GijU^ and 
u' ~ G'^^Uj, where G*-' ~ {Gij)~^ and G = dct(Gy). The six local coordinate systems 
{x^,x'^) are based on an equiangular central projection, — 7r/4 < x^,x'^ < 7r/4. The 
metric tensor for all six faces of the cube is 

1 + tan^ x^ — tan x^ tan x'^ 



G 



cos^ x^ cos^ 



tanx^ tana;^ 



(3.1) 



where ?' = (1 + tan^ + tan^ x^)^/^ and a/G = cos^ cos^ x^. 

The shallow water equations are written in curvilinear coordinates using the fol- 
lowing definitions for divergence and vorticity 

du2 dui 
dx^ dx"^ 

and replacing them in ()2.2p . In their contravariant form the equations are 

0, (3.2) 
(3.3) 

where the Einstein summation convention is used for repeated indices. Additional 
information on the equation set can also be found in [41| . 

3.1.2. Spatial and temporal discretization. The equations of motion are 
discretizcd in space using the Pjv — spectral element method as in [58]. The 

cubed-sphere is partitioned into K elements Vl^ in which the dependent and inde- 
pendent variables are approximated by tensor-product polynomial expansions. The 
velocity on element Uf^ is expanded in terms of the A^-th degree Lagrangian inter- 
polants hi 

N N 

v^(x) = EE 4 h,{v\^)) (3.4) 

and the geopotential is expanded on the same element using the (TV — 2)-th degree 
interpolants hi 

N-l N-1 

= E E hj{rf{^)) (3.5) 

i=l J=l 

where x (^'^(x), 77'^(x)) is an affine transformation from the element il'^ on the 
cubed sphere to the reference element [—1,1] x [—1,1] pictured in Fig. 13.11 A weak 
Galerkin formulation results from the integration of the equations with respect to 
test functions and the direct evaluation of inner products using Gauss-Legendre and 
Gauss-Lobatto-Lcgendre quadrature. The positions of the corresponding GL and 
GLL quadrature points within each spectral element with mapped coordinates ^ and 
T] are depicted in Fig. 13.11 The GLL points for the co-located velocity components 
are marked by the circles, the open squares point to the GL points for scalars. The 
polynomial degrees are 7 and 5, respectively. An overview of the chosen parameters 
and base resolutions for the present study is also given in Table 13.11 

continuity of the velocity is enforced at inter-element boundaries sharing 
Gauss-Lobatto-Legendre points and direct stiffness summation (DSS) is then applied 
[60l [9] . The advection operator in the momentum equation is expressed in terms of 
the relative vorticity and kinetic energy, whereas the continuity equation relies on the 
velocity form. |62] have shown that the rotational form of the advection operator is 
stable for the P^r — PAr-2 spectral element discretization. A standard Asselin- Robert 
filtered leap-frog discretization is employed for integrating the equations of motion 
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Fig. 3.1. Positions of the GLL (circles) and GL (open squares) quadrature points on the spectral 
reference element with mapped coordinates § and r). The polynomial degrees are 7 (GLL velocity 
points) and 5 (GL scalar points). 



Table 3.1 
Overview of the resolutions in SEM. 



# Elements 
per cubed face 


Total # of 
elements 


# GL points 
per element 


# GLL points 
per element 


Approximate 
resolution 


3x3 


54 


6x6 


8x8 


5° X 5° 


4x4 


96 


6x6 


8x8 


3.2° X 3.2° 


6x6 


216 


6x6 


8x8 


2.5° X 2.5° 



in time. More advanced options, that are tailored for AMR based on the ideas in 
[55] , will be described in [M]. Since a spectral basis is employed on each clement, 
the spectrum of the advective operators must be corrected using a 2/3 rule similar to 
what is employed in global pseudo spectral models. A Boyd-Vandeven filter [611 [4] 
was used during the adaptive numerical simulations that removes the last third of the 
spectrum. 

3.1.3. Non-conforming spectral elements. The positions of the collocation 
points at the boundaries of non-conforming spectral elements do not coincide and a 
procedure is required to connect neighboring elements. Several techniques are avail- 
able including mortars [521 E] and interpolations [12]. Here, interpolations are used. 
The true unknowns at a boundary belong to the master coarse element and are passed 
to the refined slave elements by the following procedure. From Eq. p.4p . at one of 
the element boundaries Tkm = f^fc H f2,„ with ^'^ = ^™ = 1 , the trace of the solution is 



N 

v{^{e,v''))t\r,^^A{^{hv''))^T.^'k,h,{rf) (3.6) 



where hj are the Lagrange intcrpolants defined at the GLL points. To interpolate 
the master solution to the slave edge, a mapping is created from the master reference 
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element relative to the slave's reference element. Let ryj^ denote this mapping, then 

V\n{v) ^ + — 1 (3.7) 

where m S [1,2'] is the slave face number and I is the number of the refinement 
level. This spatial refinement strategy is called /i-refincnicnt which is in contrast to 
increasing the order of the polynomials (p-refinement). 
Interpolations can be expressed in matrix form as 

K],y- = /i.(^™(eo) (3.8) 

where {^iji^o Gauss-Legendre-Lobatto points used in the quadrature and in 

the collocation of the dependent variables. If v*^ are the unknowns on the edge of 
the master element, then J,'„v'' represents the master element contributions passed 
to the slave elements. Assembly of the global spectral element matrix is not viable 
on today's computer architectures with distributed memory. Instead, the action of 
the assembled matrix on a vector is performed with the help of DSS. For A, a matrix 
resulting from the spectral clement discretization of a differential operator, defined 
for the true degrees of freedom and A]^, the block diagonal matrix of the individual 
local unassembled contributions of each element as if they were disjoint, the DSS for 
conforming elements is represented by 

v^Au = v^Q^AlQu (3.9) 

where Q is a Boolean rectangular matrix representing the scattering of the true de- 
grees of freedom to the unassembled blocks (one block in per element). A non- 
conforming formulation can be obtained by replacing the scattering Boolean matrices 
g with Q = JlQ 

v'^Au^v'^{Q^jI)Al{JlQ)u. (3.10) 

The block diagonal matrix consist of the identity where the interfaces between ele- 
ments match and of the interpolation matrix ()3.8|) where interfaces are non-conforming. 
To facilitate time-stepping procedures, the matrix must be lumped by summing rows 
in the mass matrix |46| . Let L{A) represent the lumping operation. Lumping the 
global matrix A is equivalent to lumping the local matrix JJ^A^Jl and it can be 
shown that 

L{A) = L{Q^JIAlJlQ) = Q^L{jIAlJl)Q. (3.11) 

3.1.4. Adaptation principle. The adaptation algorithm developed here is par- 
allel and usable on distributed memory computers. The cubed sphere is initially tiled 
with an uniform, low resolution, mesh and each element represents the root of a quad- 
tree. While a quad-tree is a natural description for 2-dimensional mesh refinement, 
the six faces of the cubed-sphere geometry necessitates a graph data structure. For 
the quad-trees, a set of simple bit manipulation subroutines are used to provide inher- 
itance information. A graph data structure which describes the connection between 
the roots and leaves of each quad-tree is also maintained. To simplify the interface 
management, neighboring elements are restricted to be at most one level of refinement 
apart. When elements are marked for refinement, a parallel procedure verifies that a 
compatible quad-tree refinement exists, with respect to all the roots and the interface 
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restriction. The communication package is properly modified to enable correct inter- 
element communications related to the DSS procedure. This was greatly simplified by 
replacing the graph-based communication package in the model by the more generic 
one developed by |60| . If a load imbalance is detected then the partitioning algorithm 
is invoked and elements are migrated to rebalance the workload on each processor. 
Apart from the Message Passing Interface (MPI) library, the graph partitioning tool 
ParMeTiS [31] and the generic DSS libraries, no other specialized high-level library 
is used. The non-adapted spectral element model also supports an efficient space- 
filling curve loadbalancing approach which we plan to use in the adaptive model 
for future simulations. 

3.2. Finite volume (FV) shallow water model. 

3.2.1. Model description. The adaptive finite volume model is built upon 
the advection algorithm by |37| and its corresponding shallow water system |38j . 
The FV shallow water model is comprised of the momentum equation and mass 
continuity equation as shown in Eqns. (|2.ip and (|2.2p . Here the fiux-form of the 
mass conservation law and the vector-invariant form of the momentum equation are 
selected. In addition, a divergence damping term is added to the momentum equation. 

The finite- volume dynamical core utilizes a flux form algorithm for the horizontal 
advection processes, which, from the physical point of view, can be considered a 
discrete representation of the conservation law in finite- volume space. However, from 
the mathematical standpoint, it can be viewed as a numerical method for solving the 
governing equations in integral form. This leads to a more natural and often more 
precise representation of the advection processes, especially in comparison to finite 
difference techniques. The transport processes, e.g. for the height h — $/g of the 
shallow water system, are modeled by fluxes into and out of the finite control-volume 
where volume-mean quantities are predicted 



Here, the overbar () indicates the volume- mean, F and G denote the time- averaged 
ID numerical fluxes in longitudinal and latitudinal direction which are computed 
with the upstream-biased and monotonic PPM scheme (see also [6j and [40]). At 
symbolizes the time step, a stands for the radius of the Earth, the indices i and 
j point to grid point locations in the longitudinal (A) and latitudinal (0) direction 
and the index n marks the discrete time level. In addition, AA^ = (A^+i — Aj_i) 
and A(j)i = {4>j+i — 4>j^\) represent the longitudinal and latitudinal grid distances 
measured in radians. 

The advection algorithm shown in Eq. (j3.12p is the fundamental building block 
of the horizontal discretization in the FV shallow water model. It is not only used to 
predict the time evolution of the mass (Eq. (|2.2[) '). but also determines the absolute 
vorticity fluxes and kinetic energy in the momentum equation. Further algorithmic 
details of the shallow water code are described in [38] who also discussed the staggered 
Arakawa C and D grid approach utilized in the model. 

3.2.2. Block-structured adaptive mesh approach. The AMR design of the 
FV shallow water model is built upon a 2D block-structured data configuration in 
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(3.12) 



spherical coordinates. The concept of the block data structure and corresponding 
resolutions is fully described in [55]. In essence, a regular latitude- longitude grid is 
split into By x non-overlapping spherical blocks that span the entire sphere. Each 
block is logically rectangular and comprises a constant number of Ny x grid points 
in latitudinal and longitudinal direction. If, for example, 6x9 grid points per block are 
selected with ByXB^ =6x8 blocks on the sphere then the configuration corresponds 
to a 5° X 5° uniform mesh resolution. Here, we use such an initial configuration for 
selected test cases in Section All blocks are self-similar and split into four in the 
event of refinement requests, thereby doubling the spatial resolution. Coarsenings 
reverse this refinement principle. Then four "children" are coalesced into a single self- 
similar parent block which reduces the grid resolution in each direction by a factor of 
2. As in SEM neighboring blocks can only differ by one refinement level. This leads 
to cascading refinement areas around the finest mesh region. In addition, the blocks 
adjacent to the poles are kept at the same refinement level due to the application of 
a fast Fourier transform (FFT) filter in longitudinal direction for nonlinear flows. 

The FV model supports static and dynamic adaptations which are managed by 
an adaptive block-structured grid library for parallel processors [331 SI] ■ This library 
utilizes a quad-tree representation of the adapted grid which is similar to the adaptive 
mesh technique in SEM. The communication on parallel computer architectures is 
performed via MPI. Our load-balancing strategy aims at assigning an equal workload 
to each processor which is equivalent to an equal number of blocks in the FV shallow 
water model. No attempt is made to keep geographical regions on the same processor. 
This can be achieved by a space-filling curve load-balancing strategy as in [5] which 
reduces the communication overhead. Such an approach is the subject of future 
research. The refinement criteria are user selected. In particular, simple geopotential 
height thresholds or vorticity-based criteria are assessed for our dynamic adaptation 
tests. For static adaptations though, we place the fine-grid nests in pre-determined 
regions like mountainous areas or geographical patches of interest. 

Each block is surrounded by ghost cell regions that share the information along 
adjacent block interfaces. This makes each block independent of its neighbors since 
the solution technique can now be individually applied to each block. The ghost cell 
information ensures that the requirements for the numerical stencils are satisfied. The 
algorithm then loops over all available blocks on the sphere before a communication 
step with ghost cell exchanges becomes necessary. The number of required ghost 
cells highly depends on the numerical scheme. Here, three ghost cells in all horizontal 
directions are needed which are kept at the same resolution as the inner domain of the 
block. The consequent interpolation and averaging techniques are further described 
in [25] and [28] . 

The time-stepping scheme is explicit and stable for zonal and meridional Courant 
numbers |CFL| < 1. This restriction arises since the semi-Lagrangian extension of 
the [37] advection algorithm is not utilized in the AMR model experiments. This 
keeps the width of the ghost cell regions small, but on the other hand requires small 
time steps if high wind speeds are present in polar regions. Then the CFL condition 
is most restrictive due to the convergence of the meridians in the latitude-longitude 
grid. Therefore, we restrict the adaptations near the pole to very few (1-2) refinement 
levels. 

The adaptive grids of both models are compared in Fig. 13.21 Here, an idealized 
mountain as indicated by the contour lines is refined at the maximum refinement level 
of 3. This corresponds to the grid resolution 0.625° x 0.625° in the finest adaptation 
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Fig. 3.2. Adaptive grids on the sphere in (a) SEM and (b) FV. The adapted elements on 
the cubed-sphere (SEM) and adapted blocks in latitude- longitude geometry (FV) refine an idealized 
mountain as indicated by the contour lines. 

region. The adaptive elements (SEM) and blocks (FV) overlay the sphere. The 
differences between the two adaptive meshes are clearly visible. The SEM model 
(Fig. 13.2b .) utilizes a non-orthogonal cubed-sphere geometry and places the refined 
spectral elements across an edge of two cubed-sphere faces. The FV model fFig. 13.2b ) 
is formulated on a latitude-longitude grid. This leads to orthogonal blocks which are 
closely spaced in polar regions. 

4. Numerical experiments. For our AMR model intercomparison, we select 
test cases with increasing complexity from the standard [63j shallow water test suite. 
The [63j test suite assesses scenarios that are mainly characterized by large-scale 
smooth flow fields. Among them are the passive advection of a cosine bell (test case 
1), a steady-state geostrophic flow at various rotation angles (test case 2), a flow over 
an idealized mountain (test case 5) and a Rossby-Haurwitz wave with wavenumber 4 
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(test case 6). Our discussion of the adaptive models is focused on these tests. Both 
static and dynamic refinement areas are addressed that highlight the strengths and 
weaknesses of the two AMR approaches. 

Some basic statistics of the adaptive simulations can be found in Table 14.11 The 
table lists not only the test case numbers with corresponding rotation angles a and 
base resolutions, but also gives information on the AMR approach, the number of 
refinement levels, the time steps as well as the number of adapted elements (SEM) 
and blocks (FV) at certain snapshots in time. Note that the time steps are not 
optimized for performance. Rather they are selected to guarantee numerical stability. 
The model results are typically evaluated via normalized h, I2 and l^c error norms. 



Table 4.1 

Statistics for the adaptive SEM and FV simulations. 



Test 




Base 
resolution 




AMR 




SEM 




FV 


At (a) 


# Elements 


At (s) 


# Ulocks 


1 




5" X 5" 


3 


dynamic 


10 


322 (day 9) 




201 (day 9) 


1 


90° 


5° X 5° 


3 


dynamic 


10 


242 (day 3) 




480 (day 3) 


2 


45° 


2.5° X 2.5° 


2 




10 


324 


200 


288 


5 




5° X 5° 


3 


dynamic 


20 


1448 (day 15) 


138 


744 (day 15) 


6 




2.5° X 2.5° 


1 




15 


312 


225 


288 


6 




3.2° X 3.2° 


1 




15 


176 







The definitions of the norms are given in |63j . For the adaptive FV shallow water 
model the computation of the norms is also explained in |28j . 

4.1. Passive advection of a cosine bell. We start our intercomparison of the 
two AMR approaches with the traditional solid body rotation of a cosine bell around 
the sphere (see test case 1 in [G^ for the initial conditions). This test evaluates the 
advective component of the numerical method in isolation. It challenges the spectral 
element method of SEM due to possible over- and undershoots of the transported 
feature. They are due to the so-called Gibbs phenomenon which is characteristic 
for spectral approaches [TH [HTJ |3] . Several values for the rotation angle a can be 
specified. It determines the angle between the axis of the solid body rotation and 
the poles in the spherical coordinate system. Here, we mainly report results for 
a = 45° which advects the cosine bell slantwise over four corner points and two edges 
of the cubed-sphere. In addition, results for a = 0° (transport along the equator) 
and a = 90° (transport over the poles) are discussed. The length of the integration 
is 12 days which corresponds to one complete revolution of the cosine bell around 
the sphere. The base resolution for both shallow water models is a 5° x 5° coarse 
mesh with additional initial refinements surrounding the cosine bell. In particular, 
the base resolution in SEM is represented by 54 elements (see Table 13. ip whereas the 
adaptations in FV start from 6x8 initial blocks. We test the refinement levels 1, 2, 3 
and 4 which are equivalent to the resolutions 2.5° x 2.5°, 1.25° x 1.25°, 0.625° x 0.625° 
and 0.3125° x 0.3125° within the refined areas. Note that these grid spacings only 
represent approximate resolutions for SEM because of the clustering of the collocation 
points near the element boundaries (see Fig. 13. ip . 

In both models we refine the grid if the geopotential height field of the cosine 
bell exceeds h > 53 m. This value corresponds to approximately 5% of the initial 
peak value with hmax = 1000 m. Of course, other refinement thresholds are also 
feasible. The threshold has been chosen since the refined area now tightly surrounds 
the cosine bell in the regular latitude-longitude grid (model FV). In SEM though, the 
adaptations are padded by a one-element wide halo. This leads to a broader refined 
area in comparison to FV but less evaluations of the refinement criterion. For the FV 
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model, the refinement criterion is evaluated at every time step. The grid is refined 
up to the maximum refinement level whenever the criterion is met. Grid coarsenings 
are invoked after the cosine bell left the refined area and consequently, the criterion 
is no longer fulfilled. The time steps in the FV simulations are variable and match 
a CFL stability criterion of |CFL| < 0.95 at the finest refinement level. This setup 
corresponds exactly to the adaptive advection experiments with FV in [28j . For the 
SEM model, the time-step is constant and restricted by the CFL number computed 
on the finest grid. 

Figure 14.11 shows snapshots of the cosine bell with rotation angle a = 45° and 
three refinement levels (0.625° x 0.625°) at day 3, 6, 9 and 12. The SEM model 
with its refined spectral elements on the cubed sphere is depicted in the left column 
(Fig. I4.1b -d). the right column (Fig. I4.1b -g) displays the FV model with its block- 
structured latitude-longitude grid. Note that each spectral element (SEM) contains 
additional 6x6 GL grid points, whereas each block (FV) consists of 6 x 9 grid points 
in the latitudinal and longitudinal direction, respectively. In both models the initial 
state (not shown) resembles the final state (day 12) very closely. It can clearly be 
seen that the adapted grids of both models reliably track the the cosine bell without 
visible distortion or noise. 

The errors can be further assessed in Fig. 14.21 which shows the time traces of 
the normalized h, I2 and loo norms of the geopotential height field for rotation angle 
a = 45°. At this rotation angle the errors in SEM are an order of magnitude smaller 
than the errors in FV. This is mainly due to the operator-splitting approach in FV 
which produces the maximum error at a = 45°. It is apparent that SEM transports 
the cosine bell rather smoothly across the corners of the cubed mesh. The errors in 
FV are considerably lower for rotation angle a = 0°. This is shown in Table l4?2l that 
documents the a = 0° normalized height errors as well as the absolute maximum 
and minimum of the height field after one full revolution for the refinement levels 1 
through 4. The errors of both models at any refinement level are now comparable to 
each other. In addition, the table shows that SEM introduces undershoots during the 
simulation which are documented by the negative minimum height values hmin ■ This 
is not the case for the monotonic and conservative FV advection algorithm. Here it 
can also be seen that the maximum height of the cosine bell hmax in FV is affected 
by the monotonicity constraint which acts as a nonlinear dissipation mechanism. It 
consequently reduces the peak amplitude of the cosine bell, especially at lower resolu- 
tions when the bell is not well-resolved. This decrease in maximum amplitude is less 
profoimd in SEM. 

Another snapshot of the cosine bell advection test at day 3 with three refinement 
levels is shown in Fig. 14.31 Here, the rotation angle is set to a = 90° which directs the 
bell straight over the poles. At day 3, the bell reaches the North Pole and both models 
refine the grid reliably over this region. Nevertheless, a distinct difference between the 
two models is the total number of refined elements (SEM) and blocks (FV) needed at 
high latitudes. While SEM (Fig. 14.3b ) keeps the total number of adapted elements 
small in the chosen cubed-sphere geometry, the FV model (Fig. I4.3b l needs to refine 
a large number of blocks due to the convergence of the meridians in the latitude- 
longitude mesh. This leads to an increase in the computational workload in the 
adaptive FV experiment due to the increased number of blocks and small time steps 
needed in polar regions. In this FV run the time step varies between approximately 
At = 600 s and At = 10 s which strongly depends on the position of the cosine bell 
P5] . The SEM model though utilizes a constant time step of At = 10 s. Further 
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90 180 270 360 90 180 270 360 
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Fig. 4.1. Snapshots of the cosine bell (test case 1) with rotation angle a = 45° and three 
refinement levels at day 3, 6, 9 and 12. The finest grid spacing is 0.625° X 0.625°, the adapted 
spectral elements (SEM) and blocks (FV) are overlaid. Left column (a-d) shows the SEM model, 
right column (e-h) depicts the FV model. The contour interval is 100 m, the zero contour is omitted. 



statistics for test case 1 are listed in Table 14.11 In particular, for rotation angle 
a — 45° the number of refined blocks in FV is smaller than the number of refined 
elements in SEM. This effect is mainly due to the additional halo region in SEM. 



4.2. Steady-state geostrophic flow. Test case 2 is a steady-state zonal geostrophic 
flow, representing a balance between the Coriolis and gcopotential gradient forces in 
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Fig. 4.2. Time traces of the normalized li, I2 and loo geopotential height error norms for the 
cosine bell advection test (test case 1) with a = 45°. The (a) SEM and (b) FV advection models 
utilized three refinement levels that correspond to the finest grid resolution 0.625° X 0.625°. Note 
that the scales in (a) and (b) differ by a factor of 10. 



the momentum equation. The initial velocity field on the sphere is given by 

u = Mo ( cos 6 cos a + cos A sin sin a ) 
V = — uo sin A sin a . 

where a is the angle between the axis of solid body rotation and the polar axis. The 
analytic geopotential field $ = gh is specified as 

$ = $0 ^ ^ afiuo + ^ ^ ( ~ A cos 9 sin a + sin 9 cos a )^ . 

a is the radius of the earth and 17 is the rotation rate. Parameter values are specified 
to be uq ~ 27ra/(12 days) and $0 = 2.94 x 10* m^/s^. These fields also represent the 
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Table 4.2 

Normalized height error norms and characteristics for the advection of a cosine bell (test case 
1) after one revolution (day 12) at different refinement levels and rotation angle cx = 0° . The first 
column indicates the resolution of the grid within the refined area. Both the SEM model (top) and 
FV model (bottom) are shown. 



Resolution 


h 


h 




h{in) max/min 


SEM 










2.5° X 2.5° 


0.0503 


0.0269 


0.0195 


991.6/-15.1 


1.25° X 1.25° 


0.0085 


0.0056 


0.0057 


997.5/-4.2 


0.625° X 0.625° 


0.0019 


0.0014 


0.0019 


999.1/-1.1 


0.3125° X 0.3125° 


0.0008 


0.0006 


0.0015 


999.7/-0.9 


FV 










2.5° X 2.5° 


0.0341 


0.0301 


0.0317 


949.1/0 


1.25° X 1.25° 


0.0097 


0.0103 


0.0150 


984.2/0 


0.625° X 0.625° 


0.0016 


0.0021 


0.0044 


995.0/0 


0.3125° X 0.3125° 


0.0003 


0.0005 


0.0014 


998.4/0 



analytic solution of the flow. The Coriolis parameter is given by 
f = 2Vl{~ cos A cos 6 sin a + sin cos a ) . 

We primarily use this large-scale flow pattern to assess the characteristics of the 
fine-coarse mesh interfaces in two statically adapted grid configurations. Both models 
SEM and FV utilize non-conforming meshes with a resolution jump by a factor of two 
in each direction at the interfaces. Note that the FV model requires interpolations in 
the ghost cell regions of neighboring blocks whenever the resolution is changed. This 
is not the case for SEM. In general, inserting a refined patch of elements (SEM) or 
blocks (FV) in a random location should result in either no changes in the error, if 
the flow is completely resolved, or in a small decrease in the overall error. Inserting 
the patch in a strategic location might even lower the error more significantly. 

Here, we start our simulations from a uniform 2.5° x 2.5° base grid which is given 
by 216 spectral elements with 6x6 GL points in SEM. In FV, this corresponds to 
the block data structure consisting of 12 x 16 = 192 blocks with 6x9 grid points in 
latitudinal and longitudinal direction, respectively. Two refinement levels are utilized 
which lead to the finest mesh spacing 0.625° x 0.625° in both models. 

In the first configuration, a refined patch of size 45° x 30° (longitudes x lati- 
tudes) is centered at (180°E, 45°N). This patch spans the domain (157.5°E,30°N)- 
(202.5°E,60°N) which is shown by the dotted contours of the adapted FV blocks in 
Fig. 14.41 Here, the geopotential height field of the FV simulation at day 14 is displayed 
which is visually indistinguishable from the initial state. In the second configuration 
a patch of identical size is centered at (135°E, 30°N). This covers the region between 
(112.5°E,15°N) - (157.5°E,45°N) which is characterized by strong gradients in the 
geopotential height field (see also the discussion in [26]). The two locations give 
insight into the dependency of the errors on the position of the refined patch. 

Both statically adapted configurations are integrated for 14 days with a = 45°. 
As mentioned before this rotation angle represents the most challenging direction for 
both models due to the choice of the cubed-sphere geometry in SEM and the operator- 
splitting approach in FV. The results for the I2 normalized geopotential height errors 
for SEM and FV arc reported in Fig. 14.51 It is expected that the errors in SEM 
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Fig. 4.3. North-polar stereographic projection of the cosine bell (test case 1) with rotation angle 
a = 90° and three refinement levels at day 3 for (a) SEM and (b) FV. The finest grid spacing is 
0.625° X 0.625°, the adapted spectral elements (SEM) and blocks (FV) are dotted. The contour 
interval is 100 m, the zero contour is omitted. 



are lower than in FV which is confirmed by Figs. I4.5b -b. This is due to the spectral 
convergence (SEM) of the smooth solution which is infinitely differentiable. In general, 
such spectral convergence can not be archived by grid point models which typically 
exhibit higher error norms. Here it is interesting to note that the two refined patches 
in SEM lead to a slight decrease in the error in comparison to the uniform-mesh run 
(Fig. 14.5k ). The errors are independent of the location of the refined area. This is in 
contrast to the error norms in FV (Fig. 14.5b ) . The refined patches introduce slight 
disturbances of the geostrophic balance and non-divcrgcnt wind field which now cause 
the error norms to increase in comparison to the uniform-mesh run. The increase is 
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90 180 270 360 

Longitude 

Fig. 4.4. Geopotential height field for test case 2 with a = 45° at day I4 as simulated with the 
FV model. The statically adapted blocks with 2 refinement levels are overlaid (dotted contour). The 
refined patch is centered at (ISCPE, 45° N), the finest grid spacing is 0.625° X 0.625°. The contour 
interval is 200 m. 

sensitive to the location of the refined pateh. In particular the errors in FV grow 
faster if the patch intersects the strong gradient regime in the geopotential height 
field (centered at (135°E, 30°N)). Despite this characteristic, the error is still in the 
expected range for grid point based models, as for example compared to [59] . The 
increase in error is mainly triggered by interpolations in the ghost cell regions (see 
also [251I2I]) and will be further investigated. 

4.3. Flow over an idealized mountain. Test case 5 of the [63] test suite is a 
zonal flow impinging on a mountain. The mean equivalent depth of the atmosphere 
is set to ho = 5960 m. The mountain height is given by hg = hs„{l — f/R), where 

= 2000 m, R = 7r/9, and = min[i?2, (A - + (6* - O^f]. The center of 
the mountain is located at Ac = 37r/2 and 9c — 7r/6 in spherical coordinates. The 
test case is integrated for 15 model days with a 5° x 5° base grid as in test case 1. 
In addition, three levels of static refinements are introduced where the height of the 
mountain is greater than m. The corresponding initial grids for SEM and FV are 
shown in Figs. [46k. [46b and[3?2l 

During the course of the simulation dynamic refinements with three refinement 
levels are applied whenever the absolute value of the relative vorticity C, is greater than 
2 X 10~^ s~^. This refinement criterion is evaluated every two hours. Grid coarsenings, 
on the other hand, are invoked in regions where the threshold is no longer met. Figure 
14 .61 shows the time evolution of the geopotential height field for both models at days 0, 
5, 10 and 15. The adapted elements on the cubed sphere (SEM) and adapted blocks 
(FV) are overlaid. It can be seen that both models track the evolving lee-side wave 
train reliably and place most of the adaptations in the Southern Hemisphere at day 
15. Some differences in the refinement regions are apparent. Overall, the refinements 
in SEM cover a broader area because a one-element wide halo was enforced around 
the regions marked for refinement. This option reduces the number of adaptation 
cycles but on the other hand increases the overall workload. The number of adapted 
elements in SEM and adapted blocks in FV at day 15 is also documented in Table [4Tl 
They almost differ by a factor of 2. Note, that SEM employs a 20 s time step, the 
time step for FV is 138 s. 

The results can be quantitatively compared via normalized error norms which 
are shown in Fig. 14.71 An analytic solution is not known. Therefore, the normalized 
error metrics are computed by comparing the SEM and FV simulations to a T213 
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Fig. 4.5. Time traces of the normalized I2 geopotential height error norms for the steady state 
geostrophic flow (test case 2) with rotation angle a = 45°. Two simulations with static refinement 
regions centered at (180° E, 45°) and (135° E, 30°) are compared to a uniform 2.5° X 2.5° resolution 
run. The adaptive (a) SEM and (b) FV model runs start with a 2.5° X 2.5° base grid with two 
refinement level (finest grid resolution is 0.625° X 0.625°. J 



spectral transform reference solution provided by the National Center for Atmospheric 
Research (NCAR) [29l|30]. The latter is available as an archived netCDF data set with 
daily snapshots of the spectral transform simulation. The T213 spectral simulation 
utilized a Gaussian grid with 320 x 640 grid points in latitudinal and longitudinal 
direction which corresponds to a grid spacing of about 55 km at the equator. Figure BTTI 
compares the normalized I2 height errors of the adaptive runs to uniform-resolution 
simulations. Both models are depicted. It is interesting to note that the errors 
in SEM (Fig. 14.7b ) already converge within the uncertainty of the NCAR reference 
solution (see [57]) at the uniform resolution 2.5° x 2.5°. The SEM adaptive run 
matches this error trace very closely. On the other hand, the errors in FV (Fig. l4.7l b) 
converge within the uncertainty of the reference solution at the finer uniform resolution 
0.625° X 0.625°. Here, the error trace of the FV adaptive simulation resembles the 
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Fig. 4.6. Snapshots of the geopotential height field (test case 5, flow over a mountain) with 
three refinement levels at day 0, 5, 10 and 15. The finest grid spacing is 0.625° X 0.625°, the adapted 
spectral elements (SEM) and blocks (FV) are overlaid. Left column (a-d) shows the SEM model, 
right column (e-h) depicts the FV model. The refinement criterion is |C| > 2 X 10^^ s^^ . The 
contour interval is 100 m. 

2.5° X 2.5° uniform run despite the higher resolution in the refined areas. It indicates 
that the coarser domains in FV stiU contribute considerably to the global I2 error 
measure. In addition, the interpolations in the ghost cell regions add to the error 
norms. 

4.4. Rossby-Haurwitz wave. The initial condition for test case 6 of the [63] 
test suite is a wavenumbcr 4 Rossby-Haurwitz wave. These waves are analytic solu- 
tions to the nonlinear non-divergent barotropic vorticity equation, but not closed-form 
solutions of the barotropic shallow water equations. However, in a shallow water sys- 
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Fig. 4.7. Time traces of the normalized I2 geopotential height error norms for the flow over a 
mountain (test case 5). The adaptive simulations with three refinement levels (0.625° X 0.625° at 
the finest level) is compared to uniform resolution runs. 



tern the wave pattern moves from west to east without change of shape during the 
course of the integration. The initial conditions arc fuUy described in |63| , the orog- 
raphy field is set to zero. The Rossby-Haurwitz wave with wavenumber 4 exhibits 
extremely strong gradients in both the geopotential and the wind fields. This test 
is especially hard for the FV adaptive grid simulations due not only to the strong 
gradients but also to the dominant 45° transport angle of the flow in midlatitudes. 
This challenges the underlying operator splitting approach in FV and accentuates 
even minor errors. 

Both shallow water models are integrated for 14 days on a 2.5° x 2.5° base grid. 
This base resolution is identical to test case 2. In addition, static refinements at 
refinement level 1 (1.25° x 1.25°) are placed within eight pre-determincd regions of 
interest. In particular, the grid is refined where the initial meridional wind field 
is u < —60 m s~^ which leads to an almost identical number of refined elements 
and blocks in the two models (Table HTT|) . This rather arbitrary refinement criterion 
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is intended to test whether the wavenumber 4 pattern moves smoothly through the 
refined grid patches. As an aside, minor improvements of the error norms are expected. 
The daily simulation results of SEM and FV are compared to the NCAR reference 
solution [29l[30]. As for test case 5 (flow over the mountain), the latter is available 
as an archived netCDF data set that contains the daily results of an NCAR spectral 
shallow water model at the resolution T213 (w 55 km). 

Figures I4.8b -b show the geopotential height field of the adaptive SEM and FV 
model runs at day 7. It can clearly be seen that both SEM and FV maintain the 
wavenumber 4 pattern of the height field rather well while moving through the refined 
patches. Here, the results can be visually compared to the NCAR reference solution 
(Fig. 14.8b ). No noise or distortions arc visible at this stage. However, FV develops 
slight asymmetries in the height field at later times. They originate at the interfaces of 
the refined areas and slightly disturb the wave field over the course of the integration. 
These perturbations can again be traced back to ghost cell interpolations in the FV 
model that now lead to a slight increase in the I2 errors norm. The effect is amplified 
by the very strong gradients in the flow field and the dominant 45° transport angle. 
The time steps for these simulations are 15 s and 225 s for SEM and FV, respectively. 

The error norms for both models are quantitatively compared in Fig. 14.91 The 
figure shows the time evolution of the normalized I2 height error norms in comparison 
to fixed-resolution SEM and FV model runs. In particular, two different base resolu- 
tions are shown for SEM. These are the regular 2.5° x 2.5° base grid (216 elements) 
and a coarser 3.2° x 3.2° base mesh (96 elements). In both SEM cases, 6x6 GL 
quadrature points per spectral element are used. The coarser simulation is added 
since SEM already converges within the uncertainty of the reference solution |57j on 
the 2.5° base grid. Figure l4!9l shows that SEM exhibits smaller errors than FV at the 
same uniform resolution. At day 12. the errors in the FV run are rather similar to 
the coarser 3.2° SEM simulation. With respect to the static adaptations. SEM does 
not show an increase in the error in comparison to FV. Instead, the SEM errors in 
the statically adaptive runs are almost identical to the uniform-mesh simulation or 
slightly diminish at the end of the forecast period. 

5. Conclusions. In this paper, two shallow water models with adaptive mesh 
refinement capabilities are compared. The models are an interpolation-based spec- 
tral element model (SEM) on a cubed-sphere grid and a conservative and monotonic 
finite volume model (FV) in latitude-longitude geometry. Both adaptive mesh ap- 
proaches utilize a quad-tree AMR technique that reduces the mesh spacings by a 
factor of two during each refinement step. Coarsenings reverse this adaptations prin- 
ciple. Then four "children/leaves" are coalesced which doubles the grid distances. In 
SEM, the refinement strategy targets the spectral elements which contain additional 
Gauss-Legendre and Gauss-Lobatto-Legendre collocation points for scalar and vector 
components. In FV, the adaptations are applied to a block-data structure in spherical 
coordinates which consists of a fixed (self-similar) block size. These blocks are sur- 
rounded by ghost cell regions which require interpolation and averaging procedures at 
fine-coarse mesh interfaces. No ghost cells areas arc needed in SEM. In both models 
neighboring elements or blocks can only differ by one refinement level. This leads to 
cascading and non-conforming refinement regions around the finest mesh. 

The models are compared via selected test cases from the standard [63j shallow 
water test suite. They include the advection of a cosine bell (test 1), a steady-state 
geostrophic flow (test 2). a flow over a mountain (test 5) and a Rossby-Haurwitz 
wave (test 6). Both static and dynamics adaptations are assessed which reveal the 
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Fig. 4.8. Snapshots of the geopotential height field of the Rossby-Haurwitz wave (test case 6) 
at day 7 for (a) SEM, (b) FV and (c) the NCAR spectral model (reference solution). SEM and FV 
utilize an adaptive grid with one level of static adaptations. The base resolution is 2.5° X 2.5° and 
1.25° X 1.25° within the refined region. The adapted elements (SEM) and blocks (FV) are overlaid. 
The contour interval is 250 m. 
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Fig. 4.9. Time traces of the normalized I2 geopotential height error norms for the Rossby- 
Haurwitz wave (test case 6). The statically adaptive simulations with one refinement levels ('1. 25° X 
1.25° at the finest level) are compared to SEM and FV uniform resolution runs. 
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strengths and weaknesses of the AMR approaches. The AMR simulations show that 
both models successfully place static and dynamic adaptations in local regions without 
requiring a fine grid in the global domain. The adaptive grids reliably track the user- 
selected features of interests without visible distortions or noise at the mesh interfaces. 
In particular, two dynamic refinement criteria were evaluated. Among them were a 
simple geopotential height threshold (test 1) as well as the magnitude of the relative 
vorticity (test 5). The latter successfully steered the SEM and FV refined grids into 
the Southern Hemisphere at the end of the 15-day forecast period. In addition, two 
static AMR configurations in user-determined regions of interest were tested. They 
confirmed that the fiows move smoothly through the refined areas in both SEM and 
FV. Nevertheless, the FV simulations showed that small errors originate at fine-coarse 
grid interfaces. These errors were due to the interpolation and averaging mechanisms 
in the ghost cell regions of the block-data structure. This was not the case in SEM 
which exhibited a decrease in error whenever adaptations were introduced. 

Overall, the number of dynamically refined elements or blocks was comparable 
in both models if the adaptations were confined to the equatorial or midlatitudi- 
nal regions. In polar regions though, the number of refined blocks in FV exceeded 
the number of refined elements in SEM considerably due to the convergence of the 
meridians in the latitude-longitude grid. In general, it was shown that SEM ex- 
hibited smaller errors than FV for all test cases at identical resolutions. This is 
expected for SEM's high-order numerical technique despite its non-conservative and 
non-monotonic nature. The latter can cause spurious oscillations and negative values 
for positive-definite fields (test 1). In contrast, the FV technique is mass-conservative 
and monotonic which, on the other hand, introduces numerical dissipation through 
its monotonicity constraint. The dissipation is reduced at high resolutions. Then the 
SEM and FV simulations resemble each other very closely. 

The problem of conservation and the removal of some of the oscillations in SEM 
can be addressed by using the discontinuous Galerkin (DG) formulation. Future 
developments will include an evaluation of the DG method as well as efficient time- 
stepping techniques [HJ [SSI [S5l \SM avoid the very small time steps for numerical 
stability. In addition, tests involving strong jets such as in [15] and [35], effects 
of orography and a full GCM comparison will be the subject of interesting future 
research. 
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